home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 7 / Apprentice-Release7.iso / Environments / MacMETH 3.2.4 / More Examples / Hennessy4.MOD < prev    next >
Text File  |  1996-06-20  |  4KB  |  213 lines
  1. MODULE Hennessy4;
  2.  
  3. FROM Storage IMPORT ALLOCATE;
  4. FROM SYSTEM IMPORT VAL, TSIZE;
  5. FROM SYSTEM IMPORT REG, SETREG;
  6. FROM InOut IMPORT WriteLn, WriteString, WriteInt, Read, OpenOutput, CloseOutput;
  7.  
  8. CONST
  9.     fftbase = 0.0 (* 1.11 *);
  10.     fpfftbase = 4.44;
  11.  
  12.     (* fft *)
  13.     fftsize = 256;
  14.     fftsize2 = 129;
  15.  
  16.  
  17. TYPE
  18.     (* FFT *)
  19.     complex = RECORD
  20.         rp, ip: REAL
  21.     END;
  22.     carray = ARRAY [0..fftsize+1] OF complex ;
  23.     c2array = ARRAY [0..fftsize2+1] OF complex ;
  24.  
  25.     Proc = PROCEDURE;
  26.  
  27. VAR
  28.     fixed,floated: REAL; ch: CHAR;
  29.  
  30.     (* FFT *)
  31.     z, w: carray;
  32.     e: c2array;
  33.     zr, zi: REAL;
  34.  
  35.     (* global *)
  36.     seed: LONGINT;
  37.  
  38.  
  39. (* global procedures *)
  40.  
  41. PROCEDURE Getclock (): LONGINT;
  42.     TYPE P = POINTER TO LONGINT;
  43.     VAR ticks: P; tk: LONGINT;
  44. BEGIN    ticks := VAL(P, 16AH);
  45.     tk := ticks^; RETURN TRUNCD(FLOATD(tk) * (1000.0D0/60.0D0) + 0.5D0)
  46. END Getclock;
  47.  
  48. PROCEDURE Initrand ();
  49. BEGIN seed := 74755D
  50. END Initrand;
  51.  
  52. PROCEDURE Rand (): LONGINT;
  53. BEGIN
  54.     seed := (seed * 1309D + 13849D) MOD 65535D;
  55.     RETURN (seed);
  56. END Rand;
  57.  
  58.  
  59.  
  60.     PROCEDURE Cos (x: REAL): REAL;
  61.     (* computes cos of x (x in radians) by an expansion *)
  62.         VAR i, factor: LONGINT;
  63.             result,power: REAL;
  64.     BEGIN
  65.         result := 1.0; factor := 1D;  power := x; i := 2;
  66.         WHILE i <= 10D DO
  67.             factor := factor * i;  power := power*x;
  68.             IF i MOD 2D = 0D THEN
  69.                 IF i MOD 4D = 0D THEN result := result + power/FLOAT(factor)
  70.                 ELSE result := result - power/FLOAT(factor)
  71.                 END
  72.             END;
  73.             INC(i)
  74.         END ;
  75.         RETURN result
  76.     END Cos;
  77.  
  78.     PROCEDURE Min0( arg1, arg2: LONGINT): LONGINT;
  79.     BEGIN
  80.         IF arg1 < arg2 THEN RETURN arg1
  81.         ELSE RETURN arg2
  82.         END
  83.     END Min0;
  84.  
  85.     PROCEDURE Uniform11(VAR iy: LONGINT; VAR yfl: REAL);
  86.     BEGIN
  87.         iy := (4855D*iy + 1731D) MOD 8192D;
  88.         yfl := FLOAT(iy)/8192.0;
  89.     END Uniform11;
  90.  
  91.     PROCEDURE Exptab(n: LONGINT; VAR e: c2array);
  92.         VAR theta, divisor: REAL; h: ARRAY [0..26] OF REAL;
  93.             i, j, k, l, m: LONGINT;
  94.     BEGIN
  95.         theta := 3.1415926536;
  96.         divisor := 4.0; i:=1D;
  97.         WHILE i <= 25D DO
  98.             h[i] := 1.0/(2.0*Cos( theta/divisor ));
  99.             divisor := divisor + divisor;
  100.             INC(i)
  101.         END;
  102.         m := n DIV 2D ;
  103.         l := m DIV 2D ;
  104.         j := 1D ;
  105.         e[1].rp := 1.0 ;
  106.         e[1].ip := 0.0;
  107.         e[l+1D].rp := 0.0;
  108.         e[l+1D].ip := 1.0 ;
  109.         e[m+1D].rp := -1.0 ;
  110.         e[m+1D].ip := 0.0 ;
  111.         REPEAT
  112.             i := l DIV 2D ;
  113.             k := i ;
  114.             REPEAT
  115.                 e[k+1D].rp := h[j]*(e[k+i+1D].rp+e[k-i+1D].rp) ;
  116.                 e[k+1D].ip := h[j]*(e[k+i+1D].ip+e[k-i+1D].ip) ;
  117.                 k := k+l ;
  118.             UNTIL ( k > m );
  119.             j := Min0( j+1D, 25);
  120.             l := i ;
  121.         UNTIL ( l <= 1D );
  122.     END Exptab;
  123.  
  124.     PROCEDURE Fft( n: LONGINT; VAR z, w: carray; VAR e: c2array; sqrinv: REAL);
  125.         VAR i, j, k, l, m, index: LONGINT; h: REAL;
  126.     BEGIN
  127.         m := n DIV 2D ;
  128.         l := 1 ;
  129.         REPEAT
  130.             k := 0D ;
  131.             j := l ;
  132.             i := 1D ;
  133.             REPEAT
  134.                 REPEAT
  135.                     w[i+k].rp := z[i].rp+z[m+i].rp ;
  136.                     w[i+k].ip := z[i].ip+z[m+i].ip ;
  137.                     h := e[k+1D].rp*(z[i].rp-z[i+m].rp);
  138.                     w[i+j].rp := h-e[k+1D].ip*(z[i].ip-z[i+m].ip) ;
  139.                     h := e[k+1D].rp*(z[i].ip-z[i+m].ip);
  140.                     w[i+j].ip := h+e[k+1D].ip*(z[i].rp-z[i+m].rp) ;
  141.                     i := i+1D ;
  142.                 UNTIL ( i > j );
  143.                 k := j ;
  144.                 j := k+l ;
  145.             UNTIL ( j > m );
  146.             (*z := w ;*) index := 1D;
  147.             REPEAT
  148.                 z[index] := w[index];
  149.                 index := index+1D;
  150.             UNTIL ( index > n );
  151.             l := l+l ;
  152.         UNTIL ( l > m );
  153.         i := 1;
  154.         WHILE i <= n DO
  155.             z[i].rp := sqrinv*z[i].rp ;
  156.             z[i].ip := -sqrinv*z[i].ip;
  157.             INC(i)
  158.         END
  159.     END Fft;
  160.  
  161. PROCEDURE Oscar ();
  162.     VAR i: LONGINT;
  163. BEGIN
  164.     Exptab(fftsize,e) ;
  165.     seed := 5767 ; i := 1D;
  166.     WHILE i <= LONG(fftsize) DO
  167.         Uniform11( seed, zr );
  168.         Uniform11( seed, zi );
  169.         z[i].rp := 20.0*zr - 10.0;
  170.         z[i].ip := 20.0*zi - 10.0;
  171.         INC(i)
  172.     END ;
  173.     i := 1;
  174.     WHILE i <= 20D DO Fft(fftsize,z,w,e,0.0625); INC(i) END
  175. END Oscar;
  176.  
  177. PROCEDURE Time(s: ARRAY OF CHAR; p: Proc; base, fbase: REAL);
  178.     VAR timer: LONGINT;
  179. BEGIN
  180.     timer := Getclock();
  181.     p;
  182.     timer := Getclock()-timer;
  183.     WriteString(s);
  184.     WriteInt(SHORT(timer), 8); WriteLn;
  185.     fixed := fixed + FLOAT(timer)*base;
  186.     floated := floated + FLOAT(timer)*fbase
  187. END Time;
  188.  
  189. PROCEDURE main2(i: INTEGER);
  190. BEGIN
  191.     fixed := 0.0;  floated := 0.0;
  192.     Time("FFT    ", Oscar, fftbase, fpfftbase);
  193. END main2;
  194.  
  195. PROCEDURE main;
  196. BEGIN
  197.     fixed := 0.0;  floated := 0.0;
  198.     Time("FFT    ", Oscar, fftbase, fpfftbase);
  199.     WriteLn;
  200.     main2(19);
  201. END main;
  202.  
  203. BEGIN    Initrand;
  204.  OpenOutput("H4.Mac");
  205.  WriteString("Hennessy4 mit MacMETH 3.2 : "); WriteLn;
  206.  WriteLn;
  207.     main;
  208.  CloseOutput;
  209.  WriteLn;
  210.  WriteString("any key to terminate. "); WriteLn;
  211.  Read(ch);
  212. END Hennessy4.
  213.